N AS A/CR - 2006-2 1 4432 



Dynamics Simulation Model for Space Tethers 

E.M. Levin, J. Pearson, and J.C. Oldson 

Star Technology and Research, Inc., Mount Pleasant, South Carolina 


June 2006 


The NASA STI Program Office... in Profile 


Since its founding, NASA has been dedicated to 
the advancement of aeronautics and space 
science. The NASA Scientific and Technical 
Information (STI) Program Office plays a key 
part in helping NASA maintain this important 
role. 

The NASA STI Program Office is operated by Langley 
Research Center, the lead center for NASA’s scientific 
and technical information. The NASA STI Program 
Office provides access to the NASA STI Database, the 
largest collection of aeronautical and space science 
STI in the world. The Program Office is also NASA’s 
institutional mechanism for disseminating the results 
of its research and development activities. These 
results are published by NASA in the NASA STI 
Report Series, which includes the following report 
types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant 
phase of research that present the results of 
NASA programs and include extensive data 
or theoretical analysis. Includes compilations 
of significant scientific and technical data 
and information deemed to be of continuing 
reference value. NASA’s counterpart of peer- 
reviewed formal professional papers but has 
less stringent limitations on manuscript length 
and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary 
or of specialized interest, e.g., quick release 
reports, working papers, and bibliographies 
that contain minimal annotation. Does not 
contain extensive analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. Collected 
papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or cosponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from 
NASA programs, projects, and mission, often 
concerned with subjects having substantial 
public interest. 

• TECHNICAL TRANSLATION. 
English-language translations of foreign 
scientific and technical material pertinent to 
NASA’s mission. 

Specialized services that complement the STI Program 
Office’s diverse offerings include creating custom 
thesauri, building customized databases, organizing 
and publishing research results. . .even providing 
videos. 

For more information about the NASA STI Program 
Office, see the following: 

• Access the NASA STI Program Home Page at 
http:! Iwww.sti.nasa.gov 

• E-mail your question via the Internet to 
help@sti.nasa.gov 

• Fax your question to the NASA Access Help 
Desk at 301-621-0134 

• Telephone the NASA Access Help Desk at 
301-621-0390 

• Write to: 

NASA Access Help Desk 

NASA Center for AeroSpace Information 

7121 Standard Drive 

Hanover, MD 21076-1320 

301-621-0390 


N AS A/CR — 2006-2 14432 



Dynamics Simulation Model for Space Tethers 

E.M. Levin, J. Pearson, and J.C. Oldson 

Star Technology and Research, Inc., Mount Pleasant, South Carolina 


National Aeronautics and 
Space Administration 

Marshall Space Flight Center • MSFC, Alabama 35812 


June 2006 


Available from: 


NASA Center for AeroSpace Information 
7121 Standard Drive 
Hanover, MD 21076-1320 
301-621-0390 


National Technical Information Service 
5285 Port Royal Road 
Springfield, VA 22161 
703-487-4650 


11 


TABLE OF CONTENTS 


I. INTRODUCTION 1 

II. DYNAMIC MODEL AND PERTURBATION RESPONSES 2 

III. DEVELOPMENT OF THE COMPUTATIONAL MODEL 5 

IV. C++ DEMONSTRATION OF MODEL PERFORMANCE 6 

V. RESULTS 7 

VI. CONCLUSIONS AND RECOMMENDATIONS 9 

APPENDIX A-MXER SIMULATION STUDY, PART I 11 

Nomenclature 12 

1. Formulation of the Problem 13 

2. Dynamic Model and Equations of Motion 13 

3. Gravitational Field Model 15 

4. Stationary Rotation 18 

5. Small Oscillations and Eigenforms 20 

6. Decomposition of Motion 25 

7. Gravitational Forces 30 

8. Simulation Approach 34 

References 35 

APPENDIX B— MXER SIMULATION STUDY, PART II 37 

Nomenclature 38 

1. Equations of Motion in Newtonian Form 39 

2. Equation of Motion in Minakov’s Form 41 

3. Boundary Conditions 43 

4. Motion of the Center of Mass 45 

5. Quasi-Static Tension 47 

6. Simulation Approach 49 

7. Resonant Excitation 52 

8. Sensitivity to Non-Gravitational Perturbations 53 

8.1 Method of Computation 53 

8.2 Aerodynamic Forces 54 

8.3 Ampere Forces 54 

iii 


TABLE OF CONTENTS (Continued) 


8.4 Solar Radiation Pressure 55 

8.5 Thermal Expansion 55 

8.6 Creep 56 

8.7 Mass Loss 57 

9. Estimation and Control Requirements 59 

10. Conclusions 60 

References 61 


IV 


CONTRACTOR REPORT 


DYNAMICS SIMULATION MODEL FOR SPACE TETHERS 
I. INTRODUCTION 


The Momentum Exchange Electrodynamic Reboost (MXER) system is a rotating tether about 
100 kilometers long in elliptical equatorial Earth orbit designed to catch payloads in LEO with its lower 
end and throw the payloads to geosynchronous transfer orbit (GTO) or to Earth escape by releasing 
the payloads at the top of the MXER rotation. To ensure successful rendezvous between the catcher at 
the MXER tip and a payload in LEO, a high-fidelity model of the system dynamics and control must 
be developed. This model must show that the MXER tether tip can be maneuvered precisely and at 
the exact time to rendezvous with and catch a payload, to within meters of positional error and within 
meters/second of velocity error, all at high tip acceleration. 

This report deals with an investigation of the modeling of the MXER dynamics in order to pro- 
vide an accurate enough representation of the MXER system and the environmental perturbations to 
ensure that payloads can be reliably caught and released on accurate trajectories. This allows an evalua- 
tion of the control requirements for MXER operations, the implications of operations in nonequatorial 
orbits, and the use of electrodynamic forces in multiple rendezvous control. We examined the accuracy 
of various models and the accuracy of environmental measurements required to predict the MXER tip 
position after one orbit, but did not deal with the requirements of electrodynamic reboosting. 


II. DYNAMIC MODEL AND PERTURBATION RESPONSES 


To develop effective mathematical approaches to the accurate simulation of the dynamics of a 
MXER tether system, we addressed two fundamental themes: 

1. Identification of the MXER system dynamics and its responses to various perturbations, the relative 
importance of perturbing factors and environmental unknowns, and the required accuracy of the 
dynamic model; and 

2. Based upon the understanding of the dynamic behavior of MXER, development and comparison of 
different computational methods to identify the most effective approach to the simulation of MXER 
dynamics. 

The theory of pendular motions formulated in the book “Dynamics of Space Tether Systems” 
(1993) by Beletsky and Levin was applied to a MXER system consisting of a number of point masses 
connected with massive tethers of fixed lengths. Equations for determining eigenforms and eigen- 
frequencies of a rapidly rotating tether system were derived, leading to the equations of the excitation of 
eigenforms under small perturbations. Eigenforms and eigenfrequencies were computed, and the motion 
of a the MXER system was developed as a combination of the orbital motion of the center of mass, 
rotation about the center of mass, and vibrations about the line of rotation. This approach gives a very 
effective and clear way to describe the dynamic responses to perturbations. 

A reference (unperturbed) motion near the equatorial plane was defined that results in an ideal 
rendezvous; the required integration accuracy was then determined for a typical rendezvous preparation 
period of a day or a few days. 

The perturbing gravitational forces on the tether system are partly caused by the fact that the 
tether length is generally not small. After considering a number of analytical approaches, an expansion 
in a form of a polynomial series was adopted to represent the gravitational perturbations on a tether sys- 
tem of finite length. 


The minimum number of terms in the gravitational field expansion needed to support high- 
precision modeling of MXER dynamics was considered. It was determined that while this number is not 
particularly large at the apogee, it increases dramatically at the perigee, with the amount of computation 
increasing proportionally to the second power of the required order of the gravitational model. A rela- 
tively simple practical formula was derived to determine the minimum required order of the gravita- 
tional field as a function of the current geocentric radius. A similar formula was derived to evaluate the 
minimum number of terms in the expansion of the gravitational perturbation force due to a finite (not 
small) spatial span of the MXER system. 
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Assuming linear superposition, the dynamic response to each perturbation was analyzed, including 
harmonics of the Earth gravitational field, interaction of the geomagnetic field with a sample current in 
the tether, aerodynamics, solar pressure, tether temperature variations, creep, and mass loss, to provide a 
richly detailed dynamic portrait of the system. 

Based on the dynamic response study, we described the impact of each perturbation on the rendez- 
vous accuracy and derive the accuracy requirements for each perturbation model. This analysis may also 
suggest ways to compensate for environmental unknowns with in-flight measurements. 

The model based on the modal decomposition allows us to calculate the position of the MXER 
tether ends to the required accuracy of about 1 meter after one orbit. The results of this task, including a 
detailed dynamic portrait of the system, its sensitivity to perturbations, and the model accuracy require- 
ments, provide the necessary basis for finding effective computational schemes for MXER. The complete 
analysis is shown in Appendix A. 
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III. DEVELOPMENT OF COMPUTATIONAL MODEL 


The MXER dynamic model consists of multiple point masses connected in succession with 
fixed-length uniform tethers. We investigated convergence, stability, error accumulation, minimum 
number of nodes or modes, ways to compute external forces, optimal integration, maximum time step 
and propagation speed. Special attention was paid to the stiffness of the system caused by a drastic dif- 
ference between the longitudinal and transverse wave velocities. We then compared the advantages and 
disadvantages of all methods and determined the most effective approach. 

One of the benefits of the suggested analytical formulation is that it allows us to operate with 
finite formulas, and to avoid integration of the gravitational forces along the tether on each simulation 
step; this leads to a faster computational scheme. The minimum number of expansion terms for the 
required accuracy was considered. 

Suitable ways of incorporating non-gravitational forces into the high-precision modal-based 
computational scheme were investigated. One of the goals was to avoid integration of the distributed 
forces along the tether on each step. It was shown that a compact set of pre-computed coefficients can be 
used to effectively account for distributed and concentrated forces with the required precision, assuming 
that the environmental models are accurate. In test runs, the computational overhead of adding the per- 
turbation forces was not significant. 

The second model used modified Minakov’s equations, which are a form of second order partial 
differential equations with respect to the tether tension. The advantage is that they allow easier separa- 
tion of fast longitudinal oscillations from the slow transverse oscillations of the tether, and therefore, 
they can be integrated with larger time steps. 

To independently verify the results of both modal decomposition method and Minakov’s method, 
a third model with lumped masses was developed. 

The complete results of this simulation are presented in the two appendices. Assuming that the 
environmental parameters are known to the required precision, the modal decomposition method proved 
to be accurate enough to provide 9-10 digit precision in prediction of the MXER tether position over the 
time of one orbit, providing accuracy on the order of 1 meter in the tip position. The results show that 
these equations can be integrated with a time step of about 3-4 seconds of orbital time, taking only a few 
seconds on a standard PC to propagate a solution over one orbit with 9-10 digit precision. 
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IV. C++ DEMONSTRATION OF MODEL PERFORMANCE 


Because of requirements for coding the equations of motion for MXER into NASA’s MXER 
Simulation Program, we provided assistance in adapting the equations for coding at the MSFC. The 
requirement was for the equations to be coded, checked, and documented as an additional part of the 
contract. Output was provided as a text file representing a solution with predefined parameters. 

We provided a computer source code from our results written as a Microsoft Visual C++ project 
for a Win32 console application demonstrating computational techniques involved in our high precision 
MXER simulation based on the modal decomposition. The source code was delivered in May of 2005, 
and is a separate item from this final report. 
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V. RESULTS 


The model based on the modal decomposition of motion was used to evaluate the effects of the 
various perturbations on the transverse tip response. The results are shown in the table. The MXER sys- 
tem is extremely sensitive to temperature variations and creep, with the responses to these perturbations 
being 1-2 orders of magnitude higher than the others. Creep is due to the tether stress, and will be rela- 
tively constant, but the temperature response will depend on the deterioration of the tether surface, cloud 
cover on Earth, orbital motion and the rotation angle of MXER. 


Table 1. Sensitivity to perturbations. 


Perturbation 

Parameter Variation 

Transverse Tip 

Response, Meters 



Aerodynamic 

5% 

2 

Ampere Forces 

1 mA 

3 

Solar Radiation Pressure 

5% 

1-4 

Thermal Expansion 

1 K 

163 

Creep 

2.4 m/day 

41 

Mass Loss 

5-150 g/day 

<0.001 


The MXER system exhibits another surprising aspect: there is a near resonance of the main trans- 
verse mode t/ 2 and the frequency of the gravity gradient variation, which cycles twice every rotation of 
the tether system about its center of mass, or 2Q. This is a long-term resonant effect, with the amplitude 
changing gradually over multiple perigee passages. The maximum tip response builds to about 120 
meters over five perigee passages, then decreases to nearly zero over the next five perigee passages. 

A typical response is shown in Figure 1 . This resonance is difficult to avoid by changing the spin 
rate, because the transverse frequencies are proportional to Q. 
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Figure 1. Transverse amplitude over twelve orbits. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The MXER dynamic system is very sensitive to perturbations, and is extremely sensitive to 
length variations. The main environmental uncertainty is temperature, which causes the largest changes 
in tip transverse position. It is therefore important that MXER have a very accurate estimation system, 
and a very accurate control system. This will require the design of estimation and control algorithms 
based on very accurate in-flight measurements of the distances between neighboring modules. This is a 
challenging task, but the distance variations carry the signatures of all the dynamic processes, and will 
form the basis for the control signals. 
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NOMENCLATURE 


A — first end-body 
B — second end-body 
C = center of mass of the tether system 
E — longitudinal stiffness of the tether 
Jc = moment of inertia of the tether system 
L — total tether length 
Lk = length of the k-th segment of the tether 
m A = mass of the first end-body 
m B — mass of the second end-body 
rftfc = embedded mass k 
R — geocentric radius 
R e = mean radius of the Earth 
R c ~ geocentric radius of the center of mass 
s = arclength along the unstretched tether 
T — tether tension 
t — time 

7 = tether elongation 
p — tether mass per unit length 
A = longitude 

T = unit vector along the tether hne 
T x — direction of the imaginary straight tether line 
Q — rotational angular rate 
( ‘ ) — differentiation with respect to time 
( 1 ) ~ differentiation with respect to the arclength 



1. FORMULATION OF THE PROBLEM 


In this study, we consider the dynamics of a spinning tether system in an 
elliptical orbit in application to the Momentum Exchange Electrodynamic Reboost 
system. Momentum exchange tether systems have been studied in a variety of ap- 
plications since Hans Moravec’s early publication [1]. It has recently been suggested 
that momentum exchange systems can be enhanced with electrodynamic reboost 
between payload transfers [2-4]. 

The Momentum Exchange Electrodynamic Reboost system (MXER) has a 
projected tether span of up to 100 km, and spins rapidly with a period of 6-7 min. 
It is placed in an orbit with a low perigee of about 400 km and a high apogee of 
about 8000 km. To capture a payload at a perigee rendezvous, within a window of a 
few seconds, the motion of the system has to be predicted with very high precision, 
having acceptable position errors on the order of 1 m. 

While this level of precision is routinely achieved today for conventional (non- 
tethered) satellites, it is much more difficult to achieve for a 100 km long flexible 
tether system. It is the goal of this study to investigate theoretical aspects of the 
dynamics and offer a practical approach to high precision dynamic modeling of a 
typical momentum exchange tether system. 


2. DYNAMIC MODEL AND EQUATIONS OF MOTION 


For the purposes of this study, we will assume that the momentum exchange 
system consists of two end-bodies A and B , modeled as point masses m A and m B , 
respectively, and a number of embedded masses m*,, k — 1, . . . , K , connected with 
tether segments of lengths Lk, as shown in Fig. 1. Mass m A can be counted as m 0 , 
and mass m B as m K +\. Segment Lk connects masses mj, and rrik-\- \ ■ 


Lq Li 


Lk 


m A m, m 2 ... m k m k+ , m B 


Fig. 1. Structure of the momentum exchange tether system. 
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All point masses and the masses of the tether segments are assumed to be 
constant. Tether mass per unit length can vary along the tether. 

Positions of the tether elements with respect to a non-rotating geocentric 
reference frame OXY Z are defined by the geocentric radius R as a function of the 
arclength s measured along the unstreched tether from A to B, and time t, 

R = R(s,t), 

Positions of the end masses and embedded masses are 

R-A = Rfl — R(SS;t)> Rfc = R(sfc,t). 


k 


A 



Fig. 2. Positions of the tether system elements. 


The tension vector of a perfectly flexible tether is tangent to the tether line 

T = — R\ 7 = l R 'l- (1) 

7 

where prime denotes differentiation with respect to the arclength s, and 7 is the 
local tether elongation. 

The tether tension can be expressed as a function of the elongation 7, elon- 
gation rate 7, temperature 0, and other factors, such as creep history, 

T = T(s, t, 7, 7, 0 ,...). ( 2 ) 
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Equations of motion of the tether system constitute a mix of ordinary and 
partial differential equations [5]. The motion of the end masses and embedded 
masses is described by the ordinary differential equations 

m A R^ = T a + m A g A + F a 

m B Rb — — Tb + m B g B + Fg (3) 

m k Rjt = Ta. + - T fe _ +m k g k + F k 

where dots denote differentiation with respect to time, g^, g B , and g k are the 
gravity accelerations at points A, B, and k , respectively, F A , F B , and F*. are non- 
gravitational forces acting on the end masses and embedded masses. The tether 
tension vectors are taken at the following points: T A at point A, T B at point f?, 
Tfc_ at point k of segment L k ~\ , and T*. + at point k of segment L k . 

The motion of the tether is described by the partial differential equation 

pR = T'+pg + F (4) 

where dots denote differentiation with respect to time t, and primes denote differ- 
entiation with respect to the arclength s, p is the tether mass per unit length, and 
T is the tether tension. 


3. GRAVITATIONAL FIELD MODEL 


The gravitational field is a sum of the gravitational field of the Earth and 
other celestial bodies, primarily, the Moon and the Sun. 

The geopotential U is usually represented as 


U 


I^E 

~R 


EE 

n=0 m— 0 


Re 

Ji 


{Cnm cos m\ + S nm sin mA) P nm (cos a), (5) 


where is the gravitational constant of the Earth, R E is the equatorial radius of 
the Earth, R is the geocentric radius, A is the geographical longitude eastward from 
Greenwich, a is the geocentric colatitude (the angle between the rotation axis of 
the Earth and the geocentric radius of a point, a = 0 at the North Pole), P nm are 
normalized Legendre functions, 

Pnm(cosa) = x™P™(cosa), 
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is a norm, 


= \/2n + 1, 


yt™ 1 — 


1 2(2 n + 1) {n — m)\ 
(n + m)\ 


with m = 0, 
with m > 0, 


and P™ are associated Legendre functions of degree n and order m, 


P™(*) = ( l-* 2 ) 


_ 2 \m/2 


d r 


dx m 


Pn{x ), 


P n (x) = 


d n 


2 n n! dx n 


(x> - 1)", 


P° _ 1 

-r Q — i, 

P^cosa) = cos a, 

P* (cos a) = sin a, 

^2° (cos a) = -(3 cos 2a + 1) = - (3 cos 2 a - 1), 

3 

Pi (cos a) = - sin 2a = 3 sin a cos a, 

Z 

3 

P|(cosa) = -(1 — cos 2a) = 3sin 2 a, 

Z 

P“ (cos a) = -(5 cos 3a + 3 cos a) = -(5 cos 3 a — 3 cos a), 
8 2 

3 3 

P 3 1 (cosa) = -(sin a + 5 sin 3a) = - sina(5cos 2 a — 1), 

8 2 

7~>2 / \ 1 5 , . o 

P 3 (cos a) = —(cos a — cos 3a) = 15 sin a cos a, 

15 

P 3 (cos a) = — (3 sin a — sin 3a) = 15 sin 3 a, 


The components of the gravitational acceleration g = VU are 

_ d u _ i du i du 

9R ~ dR' 9a ~ Rda' 9X ~ Psina ’ 

where <7 fl is pointing along the geocentric radius, g a is pointing southward along 
the meridian, and g\ is pointing eastward along the parallel. In the geocentric axes 
OXYZ, 

g x — ( gR cos <p + g a sin <p) cos A — g\ sin A 
9y = ( 9r cos + g a sin <p) sin A + g x cos A 
9z = 9r sin ~ g a cos <p 

where p> = tt/2 — a is the latitude. 
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To properly evaluate convergence, we should use a geopotential model of the 
maximum available degree and order. An obvious choice is the Earth Gravitv Model 
EGM96 [6], of degree and order 360. 7 

According to this model, fi E = 398600.4415 km 3 /sec 2 , R E = 6378.1363 km, 
and the normalized coefficients C nm and S nm are as follows: 


C'o.o = 1, 

C'i.o = 0, 

^ 1,1 = 0, 

C 2 ,o = -0.484165371736 • 10~ 3 , 
C 2 ,i = -0.186987635955 • 10“ 9 , 
C 2 ,2 = 0.243914352398 • 10~ 5 , 
C 3 , o = 0.957254173792 • 10 -6 , 
C 3 , i = 0.202998882184 • 10" 5 , 
C 3 , 2 = 0.904627768605 • 10~ 6 , 
C 3) 3 = 0.721072657057 • 10" 6 , 


50.0 — not used 
-S'ljO — not used 
§1,1 = 0 , 

S 2t o — not used 

5 2 .1 = 0.119528012031 • 10 -8 , 
S 2l2 = - 0.140016683654 • 10 ~ 5 , 
S 3< o — not used 

§z , i = 0.248513158716 • 10 -6 , 

5 3 . 2 = - 0.619025944205 • 10 “ 6 , 

5 3 . 3 = 0.141435626958 • 10 ~ 5 , 


<^ 360,359 = 0.183971631467 • 10 10 , S 36 0,359 = -0.310123632209 • 10 _1 °, 

^360,360 = -0.447516389678 • 10 -24 , 5 36 o,36o = -0.830224945525 • 10 -1 ° 

Gravitational perturbations caused by the Earth tides can be described by 
time dependent variations in the geopotential expansion coefficients [7], 

Gravitational perturbations caused by a celestial body of mass M p located at 
R p with respect to the geocentric axes OXY Z are described as 


gp — GM p 


/ Rp — R 

VIRp-RI* 



( 6 ) 


where G is the universal gravity constant, and R v = |R p |. 

For most practical purposes and for all celestial bodies, except the Moon, the 
linear approximation 

GM P . , 

gp — [3e p (ep,R) — R], (7) 

where e p = R P /R P is a unit vector of the direction to the celestial body, should be 
adequate [7]. 
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4. STATIONARY ROTATION 


If we disregard non-gravitational forces F in equations ( 3 )-( 4 ) and assume 
that the tether is inextensible (7 = |R'| = 1 ) and gravitational acceleration does 
not vary along the tether (g = g c ), then equations ( 3 )-( 4 ) will have a stationary 
solution 

T = T(s)t(2), R — Rc(0 + (<$ “ 6 C ) t(<), (8) 

where R c(t) and t(£) are determined from 


Rc = gc, T = —Q 2 T. 


(9) 


Here, R c is the radius- vector of the center of mass, moving as a point mass, g c is 
the gravity acceleration at the center of mass, T is a unit vector along the tether line 
rotating at a constant angular velocity Q in a fixed plane normal to O, Q = |Q|, 
s c is the arclength corresponding to the center of mass 


s c 


— — (^rnA s A + r fnB s B + ^ ] mk s k + J p s ds^j , 


( 10 ) 


and M is the total mass of the tether system. 
In the stationary motion ( 8 )— (9) , we have 


R = R c + (s - s c ) T = gc - (-s - s c ) f? 2 T. 


After substituting this relation into equations (3)-(4), we arrive at a boundary 
problem for the equilibrium tension, 

T' = —pfi 2 [s — Sc), 

T a = —Tn A 0 2 (sA — Sc), 

Tb = ttib&{sb Sc), 

Tk+ = Tk- — mfcf? 2 (sfc — s c ), 

where the last equation defines tension increments associated with the embedded 
masses m^. 

The equilibrium tension is obtained by solving the first equation of (11) with 
the initial condition T = T A given by the second equation and incrementing the 
tension at the points Sk, as described by the last equation of (11). The boundary 
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condition T — T B given by third equation is always satisfied because of the definition 
of the value s c (10) and the structure of equations (11). 

In general, the equilibrium tension increases from the value T A at the end A 
to its maximum value Tc at the center of mass C , and then drops to the value Tg 
at the other end, as shown in Fig. 3. 

To minimize mass, the tether should be tapered, with the linear density vary- 
ing along the tether, p — p(s). Ideally, the tether should be equally stressed/strained 
at all points, 


T 

E 


= 8 




( 12 ) 


where E is the longitudinal tether stiffness, and is the maximum allowed strain. 
Condition (12) can be rewritten as 


T = 


P V E ? 



(13) 


where v E is the longitudinal wave velocity in the tether. Now, we can express the 
linear density as a function of tension, p = Tv^ 2 S * 1 , and substitute this relation 
into the first equation of the boundary problem (11) to produce a boundary problem 
for an ideal equally stressed/strained tether. 


t (km 



Fig. 3. Tension profile of a tether system with (a) and without (b) payload. 




Ill practice, the tether can be made of a number of uniform segments, approx- 
imating the desired taper profile. 

Typical tension profiles are shown in Fig. 3 for a 90 km tether system with 
m A - 250 kg, m B - 11000 kg, and 8 embedded masses m k = 200 kg placed 10 km 
apart. The tethers are made of Zylon and are uniform on each segment with a 
maximum strain of S„ = 0.01. The system spins at 12 = 0.8 deg/sec. In case (a), 
the end mass A carries a 2500 kg payload, and there is no payload in case (b). 


5. SMALL OSCILLATIONS AND EIGENFORMS 


To study small oscillations of the tether system about the stationary rotation, 
we introduce a rotating reference frame C£r](, with the origin at the center of mass, 
C£ aligned with the tether line, and axis G aligned with the angular rate vector 
n, so that 

t= (1,0,0), Q = (0, 0, Q). 

With respect to the rotating axes C£rj(, the stationary motion (8)-(9) is viewed as 
a relative equilibrium. 

Equations of small oscillations of the tether system about the relative equi- 
librium are derived from (3)-(4), 

pD(8 r) = 6T', 
m A D(8r A ) = 6T a , 
m B D(8 t b ) = -8T b , 
m k D(8v k ) = 6T k+ - 

where 8r and 8 T are deviations from the relative equilibrium, and D denotes the 
linear expression 

D(8r) = 8r + 2D. x 8r + Q x (D x 8r), 
in which derivatives are calculated with respect to the rotating axes C£r](. 

Small oscillations (14) have two independent components. One is normal to 
the tether line in the plane of the stationary rotation, 

£r = (0, 77, 0), 8T = (0, Tr, ', 0), 

and the other one is normal to the rotation plane 

8v = (0, 0, C), ST = (0, 0, T('). 


There is no longitudinal component in the linear approximation for an inextensible 
tether. 



Equations of small in-plane oscillations take the form 


p{r,-Q 2 r ] )^{T7 1 , )\ 
m A (rj A - f2 2 r] A ) = (Tt]') A , 
m B (vb - (2 2 Vb) = -(Trj') B , 
m k (rjk ~ fi 2 Vk) = (Trj') k + - (Tr]') k -, 


( 15 ) 


After substituting rj = U n (s) cos(Q n t) into (15), we arrive at the following eigen 
value problem, 

(Tu'j = - P (nl + tf)u n , 

{Tu' n ) A = -m A {n 2 n + n 2 )u nA , 

(TU’ n )a = m B + 12 2 ) U nB , 

(TV n ) k+ = ( TU’ n ) k - - m k (ill + Q 2 ) U nk , 

Small out-of-plane oscillations are described by 


(16) 


p< = (T(’y, 

m A Ca = (T (') A , 

t^bCb = —{T(')b, 
m k ( k =(T(') k+ -(T(') k -. 


(17) 


After substituting ( = U n (s ) cos(fl n t) into (17), we derive the following eigenvalue 
problem, 

(tu'J = - P nlu n , 

(TUUa = -m a alu nA , 

(TU’ n ) B = m B alU nB , 

{TU' n ) k+ = (TU' n ) k --m k nl U nk . 

The eigenvalue problems (16) and (18) are very similar. They have the same 
eigenfunctions U n (s ), and their eigenfrequencies are bound by a simple relation, 


n 2 n< = nL + n 2 . 


717 ? 


(19) 


Analyzing the boundary problem (11) for the equilibrium tension and the 
eigenvalue problems (16), (18), we note that the eigenforms do not depend on the 
angular rate !?. Indeed, we can introduce a normalized equilibrium tension 


p M = 


t W 
n- 


( 20 ) 


21 



and rewrite equations (11) as 


P' = -p(s - s c ), 

Pa — ~ rr >-A{ s A — •Sc)) 

P B = m B (s B - s c ), 

Pk+ = Pk- ~ rn k (s k ~ s c ), 

while equations (16) and (18) can be reduced to 

(PU'J = -pftU n , 

( pu L)a = -m A plU nA , 

{PU' n ) B = rn B (3l U nB , 

( PU'n)k+ = ( PU'Jk - - m k Pi Unk, 


(21) 


( 22 ) 


where P n are dimensionless eigenvalues. The resulting system of equations (21)- 
(22) and its solutions do not depend on the angular rate Q. 

The eigenfrequencies Q nv and Q n ^ can be expressed through the dimensionless 
eigenvalues ft n as 


n nv = y/Pi-i n, n n( =p n n. ( 23 ) 

The eigenvalue problem (21)-(22) has two trivial solutions, 

Po = 0, U 0 = 1, (24) 

and 

P\ = 1) Ui = s - s c . (25) 

The first solution simply reduces the right and the left parts of (22) to zero, while 
the second solution reduces (22) to (21). Dynamically, solution (24) corresponds to 
perturbations of the motion of the center of mass, while solution (25) corresponds 
to variations of the angular rate and direction of the axis of the stationary rotation. 

The trivial solutions are followed by an infinite series of solutions {P n , U n }, 
n = 2, 3, ... , which constitute an orthogonal basis with the orthogonality condition 

(Ui, Uj) = 0, i j, (26) 


where 


(Vi, Vi) 


— m A(UiUj)A + TnsiUiU j)s + 


J2 m ^iUj) k 

k 



p(UiUj) ds. 


(27) 


22 





23 










( 32 ) 


Also, the following relation for the eigenvalues can be derived from (22) 

0l = ^~J*P(K) 2 ds = ci. 

Some eigenforms are shown in Fig. 4-5 for the same parameters as in case (b) 
of Fig. 3, namely, L = 90 km, m A = 250 kg, m B = 11000 kg, m k = 200 kg, k = 
1,. . . ,8, connected with 10 km long uniform tethers of the following masses: 318, 
365, 406, 438, 461, 473, 475, 473, 460 kg. The corresponding nontrivial eigenvalues 
(re = 2, ... , 19) are given below 

P 2 = 2.136, /3 a = 10.497, /? 14 = 20.182, 

Pz = 3.339, p 9 = 11.355, p 1B = 20.573, 

Pa = 4.590, p 10 = 12.414, p 16 = 22.015, 

Ps = 5.912, p^ = 15.549, p 17 = 23.714, 

Pe = 7.328, p 12 = 18.020, p ls = 25.550, 

P 7 = 8.853, p 13 = 19.320, p 19 = 27.529. 

Generally, the n-th form has re nodes. Note that while the lower modes (Fig. 4) 
show smooth curves of average collective behavior, the higher modes (Fig. 5) exhibit 
a more peculiar interplay between the tether and embedded mass dynamics. 

6. DECOMPOSITION OF MOTION 


Using eigenforms obtained from the solution of the boundary problem (21) 
and the eigenvalue problem (22) for a simplified case of a stationary rotation, we can 
represent solutions of the general equations (3)-(4) in the form of a series expansion 

OO 

R ( s > <) = X! **"(*) U M> (33) 

n—0 

where the term re = 0 (trivial form U 0 = 1) represents, in essence, the orbital motion 
of the center of mass, q 0 = Rc(t), and the term re = 1 (trivial form U\ = s — s c ) 
describes a quasi-rigid rotation of an imaginary straight tether, q x = T 1 (t), where 
is pointing along the imaginary straight tether line, so that 

OO 

R(s, t) = Rc(0 + (5 — s c ) Tj (t) + ^2 q n(t) U n (s). (34) 

71=2 
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The higher order terms n > 2 represent tether oscillations about the imaginary 
straight line drawn through the center of mass C along vector t x . Kinematically, 
they do not contribute to the displacement of the center of mass and the quasi-rigid 
rotation because of the orthogonality conditions (26). 

Note that the term with T x also describes a uniform elongation of the tether 
system. If we disregard the higher order terms n > 2, then the tether elongation 
7 = |R'| = | T i| = 7 i will be the same for all tether elements. 

Expression (33) yields similar kinematic and dynamic relations, 

00 oo 

R( 5 , t) = J2 qn(<) U n (s), R (M) = Y, q n(<) U n ( S ), (35) 

n=0 n=0 


or in terms of (34), 


R(s, t) — R c (t) + (-5 s c ) Ti(t) + 'Y, Qn{t)U n (s), 

71=2 

oo 

R(a, i) = R c (i) -f (if — Sc)Ti(() + <in(t)U n (s), 


(36) 


71=2 


where the derivatives are calculated with respect to the non-rotating geocentric 
frame OXYZ. 

We can now substitute the second equation of (35) into the general equations 
(3) (4), multiply by U n , n — 0, 1, 2, . . . , and sum over all tether elements and point 
masses. Applying the orthogonality conditions (26), we find that 


Qn, — 


Wn\ 


(Qra 4" + ®n)) 


n = 0,1,2,..., 


(37) 


ds 


where Q n is reduced to an integral of tension, 

f B 

Qn = T A U nA -T B U nB + Y( T k+-T k -)U nk + / T'Z7 n 

k 

= - f B T V n ds, 

J A 

G n is expressed through the inner product (27), 

G» = (g, U n ) 

f B 

= m A (gU n ) A + m B (gU n ) B + Y m k(gU n ) k + / p(gU n )ds , 

L j A 


(38) 


(39) 



and <D n is a generalized sum of non-gravitational forces, 

f B 

^aUtia + FbUub + ^ ] Fj, U n k + / F U n ds. (40) 

k J A 


The system of ordinary differential equations (37) is equivalent to the original 
system of mixed partial and ordinary differential equations ( 3 )-( 4 ), because it is 
derived through the general transformation of variables ( 33 ) without any simplifying 
assumptions. 

Using notations (34), and taking into account that ||f7 0 || = M and Q 0 vanishes 
because U 0 = 0 , the motion of the center of mass is described by 

R-c = — (G 0 + <D 0 ). (41) 


Taking into account that ||17i|| — J c and U[ — 1 , the quasi-rigid rotation and 
uniform elongation are described by 


Ti = 


(Qi + Gi + ®i), 


Q 


i — 


/: 


T ds. 


According to (34), the tangent to the tether line is defined as 

OO 

R' = t 1 + ^t, Ar= '£q n U' n . 

71=2 

Assuming that |utJ <£ T, , , the elongation can be approximated as 


(42) 


(43) 


T= l R 'l = V / Ki| 2 + 2(T 1 ,aT) + |ziT| 2 « 7l +(e l ,4T) = 7 1 + ^T qen U' n , (44) 

71=2 

where 

Tl K|, 7 9en, = (®lj*fn)' (4b) 

The tether tension depends not only on the elongation, but also on temper- 
ature, internal friction, creep history, and other factors. The following expression 
can be adopted within a certain range of conditions, 


T = T 0 + E(j - 7 *) + Dj, (46) 

where T 0 is a reference tension, E is the tether longitudinal stiffness, D is the 
damping coefficient, 

7 * = l0 + e(t-t 0 ) + a(G-e 0 ), ( 47 ) 
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7 o is the elongation under the reference tension T 0 at time t 0 and a reference tem- 
perature 6 > 0 , e is the creep rate, a is the thermal expansion coefficient, and 0 is the 
current temperature of the tether element. Other terms can be added as needed. 
Substituting (44) into (46), we derive 

oo 

T — Ti q en + D q en ) U n , T t = T 0 + E (ji — 7 *) + (48) 

n~ 2 

and then calculate 


T 

7 




9en + D q en 


UL 


(49) 


Retaining only linear terms in q n , the integral of tension Q„ can now be 
represented as 


Qr 


T T / 00 

= ~Ja 7 R,U n is = -J A “( T > + X>tfi)tC<fa 

k~2 

fB 00 

* -J ( t ^ + T, c ‘‘U' kjKds, 


(50) 


where ej — T x / 71 is a unit vector along the imaginary straight tether line and 


T 

c n = Ee^qen + De 1 q en + ~ (q„ - e l9en ), 

7i 


?en = (e 1 ,q n ). (51) 


Expression (50) can be significantly simplified in the ideal case of a perfectly 
tapered tether. Let us define the tether reference state T 0 (s) as a stationary rotation 
at an angular rate f2 0 free of external torques, as described in Section 4 , with the 
exception that the tether is now elastic. Let us assume that the tether is perfectly 
tapered so that all tether elements have the same elongation 7 0 under the stationary 
tension T 0 at a given temperature 0 O . Then, similarly to ( 11 ), the equilibrium 
tension profile can be determined from the boundary problem 

T' = -pfi 2 (s - s c ) 70 , 

Ta = —m A Q (5^4 — Sc) Jo, 

T B = m B fi 2 (s B - s c ) 70 , ^ 

Tk-\- Tfc— (<S/fc Sc) 7o , 

and similarly to ( 20 ), the equilibrium tension can be expressed as 

To(s) = f2 0 7o P($), (53) 
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where P(s) is the solution of the boundary problem ( 21 ). 

Within the linear elasticity region, the longitudinal stiffness of the tether can 
be related to its tension as 

B(s) = ^\=k. flj P(s), = * (54) 

With a strain on the order of 7o - 1 ~ 0.01, the coefficient k E ~ 100. 

The internal damping coefficient is usually considered directly proportional to 
the longitudinal stiffness and inversely proportional to the frequency of oscillations 
f? n , which gives 


D(s,f2 n ) 


n n 


v n n 0 p(s), 


Vn 


k E rj 

^T’ 


(55) 


where tj is a loss factor, on the order of 0.1 for braided tethers. 

After these substitutions, the first term under the integral in formula (50) 
takes the form 


T iei = QlP{s )[ l0 + k E { lt -7*) + i/i7i] ei = Q 2 0 P{s)u 1 x u 

where 

— [Tq+^Ct 1 -TO + ^iTi], (56) 

and the expansion coefficients become 

Cn = f2 0 P(s) [& E e 1 q en + u n exg en + (q n — e^n)] . 

If we assume that 7 * defined by (47) does not vary along the tether, i.e., the 
creep rate, thermal expansion coefficient, and the temperature are the same for all 
tether elements, then we can apply the orthogonality condition (28) along with the 
eigenvalue relation (32) to derive 

Qi = -Jo^u 1 x 1 , ( 57 ) 

for n = 1 and 

Q™ ~ — ll^nll [kE Gi9en + V n e 1 q en + U x (q n — e x ^erj.)] (58) 

for n — 2 , 3 ,... 

Now, we can rewrite the equation of quasi-rigid rotation (42) as 



and the equations of tether oscillations (37), n = 2 , 3 ,..., 


as 


q n +/?n^O [ fc B e l?en + ^nej^en +«l (q„ - e^en)] = -J_ (Q n + <D n ). (60) 


7 . GRAVITATIONAL FORCES 


The gravitational terms in the equations of motion (41), (59) and (60) are de- 
fined by (39). The variation of the main component of the gravitational acceleration 
g along the tether can be represented in the linear approximation in displacement 
from the center of mass as 


g = 


He R 



He 


[3e* (e B ,r) - r], 


(61) 


where He is the gravity constant of the Earth, e R = R C /R C is a unit vector along 
the geocentric radius of the center of mass R c , and 


r( 5 , t) = (s- s c ) t x ( t ) + q ra (t) U n (s ) (62) 

n = 2 

is the position vector relative to the center of mass. 

Substituting expressions (61)-(62) into (39), and applying the orthogonality 
conditions (27) and the norm definition (29), we obtain 


G 0 

G x 


- M 



J o [^(en.T,) — T t ], 

\\Vn\\^ [3 ejt ( eit ,q„)-q„], 


n ~~ 2,3,..., 


(63) 


where M is the total mass of the system (30), and J c is the moment of inertia of 
an imaginary system with an unstretched straight tether about the center of mass 
defined by (31). 

Note that in the linear approximation, the higher modes of oscillations do 
not affect the motion of the center of mass and the quasi-rigid rotation about the 
center of mass, and there is no practical need to include higher order terms in q n as 
long as typical amplitudes of tether oscillations in a momentum exchange system 
are relatively small. However, to achieve the required precision, it is necessary to 
include higher order terms in (s — s^)^. 
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The variation of the main (Newtonian) term of the gravitational acceleration 
along the line ( s - s c ) x 1 drawn through the center of mass is given by 

g — ^ E ^ _ _l l E [Rc + (s — ^c)^l] 

R 3 |Rc + (5 — S c ) Ti | 3 

where 

* = -(e«,e 1 ), e x = — , 

h 

As known in the theory of Legendre functions, 


(1 - 2are + e 2 ) 1 / 2 X P n(x) 

7 n-0 

where P n (x) are Legendre polynomials. Differentiating this equation with respect 
to x and dividing by e, we find that 


— _ f J ' E e n T ee i , 

R 2 C (1 - 2xe + e 2 ) 3 / 2 ^ 


( s ~ fo) 7i 


(l - 2x£ + £ 2 ) 3 / 2 X^+lC*)^ 

7 n=0 


where -P n are the first derivatives of the Legendre polynomials. 
After substituting this expression in (64), we obtain 


g ~~W X i e * P n+i( x ) + e i K(x)] 


Rl 




n=0 


(* - g c) 7i 


(65) 


This series expansion converges rapidly because ( s-s c )/R c is typically small, 
on the order of 0.01. To determine a minimum number of terms required for 
practical computation, consider a simple inequality 


n 



or 


> In (n/S) 

~ In {Rc/LY 


(66) 


where the left part estimates the contribution of the n-th term, L is the total tether 
length, and 6 is a maximum acceptable relative error. Generally, all terms N, 
where N is the minimum number satisfying (66), can be dropped. For a typical 
momentum exchange system, it is sufficient to retain only 5 terms at the perigee, 
and 4 terms at the apogee to achieve 10-digit precision ( 8 lO” 10 ). 

Expansion (65) can be also presented as 


g = Xsc(*-sc) n , 

n— 0 


K » = if!Zs 

c n\ d s n 


( 67 ) 
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where the derivatives with respect to s are calculated along the line (s - 3c ) Tl 
drawn through the center of mass. For the Newtonian term, according to (65), 

= [ e *-Pn+l(*) + ei P' n {x)\ (68) 

Note that expansion (67) can be treated more broadly to include all terms of 
the standard representation of the Earth gravitational field, as well as the gravita- 
tional fields of other celestial bodies. 

Using (67), the generalized gravitational forces G n (39) can be represented as 


— (( s ~ s c) k , U n ). (69) 


Applying the orthogonality conditions (27), norm definition (29), and trivial 
form properties (24)-(25), (30)-(31), we find that 



where I n are the high order moments, 


/„* = h, 

If =Ik + l, 


(70) 


In 


m A(sA s c ) n + m B (s B — 5 C ) n -(- 'y ^ mfc (sfc 

k 



s c ) n ds. (71) 


We note also that 


r 

n 


In = 0 , 


n = 2,3,. 


(72) 


The practical benefit of using expansion (69) is that the quantities can be 
precomputed to avoid integration of the gravitational force along the tether on each 
step of solution propagation. 

Another practical question is to estimate how many terms of the gravitational 
field expansion (5) must be retained to achieve the required accuracy. The problem 
is in the slow convergence of the series expansion at low altitudes when the ratio 
R E / R is close to 1. 


On average, according to the Kaula rule [7], the normalized geopotential co- 
efficients C nrn and S nrn decrease in inverse proportion to the second power of their 
degree, 


Cnm ? $r 


10 - 


n * 


(73) 


32 





o 1 2 t (h) 3 


Fig. 6. Minimum degree of the geopotential expansion for MXER. 


Taking this empirical rule into account, the following inequality roughly esti- 
mates the contribution of the higher degree terms to the acceleration of the center 
of mass of a momentum exchange tether system, 


10“ 5 

n 



^ In (nS • 10 s ) 

U ~ ]n(R c /R E ) ’ 


(74) 


Numerical experiments show that a value of 6 ~ 10 -9 can be assumed for 
MXER to achieve better than 1 m prediction accuracy over one orbit. The minimum 
degree N that satisfies (74) is shown in Fig. 6 as a function of the geocentric radius 
and flight time. We see that a lower degree model is adequate for the most part of 
the orbit, but a model of a much higher degree is needed during perigee passages. 


33 





8. SIMULATION APPROACH 


Now, we can gather all equations, and apply formulas for the gravitational 
forces derived m the previous section. The equation of motion of the center of : 

(41) takes the form 


mass 


®o 

Af ’ 


k- 2 

the equation of the quasi-rigid rotation (59) can be rewritten 

1 °° 

+ fil Ul x, = gl + -1 + p. 


(75) 


as 


C 1^2 J ° 


(76) 


and the equations of tether oscillation (37), n = 2, 3, . . . , can be reduced to 

qn+ [k E ®1 9en + (q n — e 1 g' en )] = 

g(, 


(77) 


k~ 2 


As noted earlier, very few terms g£ are needed to achieve the desired precision, 
and quantities can be precomputed. 

Equations (75)-(77) cleanly separate three dynamicahy different kinds of mo- 
tions of the momentum exchange tether system. They are compact and flexible 
in sense that the number of the gravitational terms and the number of the tether 
oscillation modes can be selected depending on the desired accuracy. This opens a 
way to building a very computationally effective simulator. 

Preliminary testing for a typical momentum exchange system showed that 
equations (75)-(77) can be integrated with a time step of about 4 seconds of orbital 
time, and it takes only a few seconds on a standard issue PC to propagate a solution 
over one orbit with 9-10 digit precision. More studies are needed to understand 
possible limitations of this approach. 
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NOMENCLATURE 


A = first end-body 
B = second end-body 
C — center of mass of the tether system 
E — longitudinal stiffness of the tether 
F = non-gravitational forces 
g = gravitational acceleration 
Jo — moment of inertia of the tether system 
L = total tether length 

Lk — length of the &-th segment of the tether 
m a = mass of the first end- body 
= mass of the second end-body 
rrik = embedded mass k 
M — total mass of the tether system 
R = geocentric radius-vector 
Re — mean radius of the Earth 
Rc = geocentric radius of the center of mass 
$ — arclength along the unstretched tether 
T — tether tension 
T = tether tension vector 
i = time 

7 = tether elongation 
p = tether mass per unit length 
T = unit vector along the tether line 
T( = direction of the imaginary straight tether line 
1? = rotational angular rate 
( ’ ) — differentiation with respect to time 
( 1 ) — differentiation with respect to the arclength 
(a,b) — scalar product of vectors a and b 
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1 . EQUATIONS OF MOTION IN NEWTONIAN FORM 


In the first part of this study [1], we considered the dynamics of a momentum 
exchange system consisting of two end-bodies and a number of power stations con- 
nected with tether segments. The end-bodies A and B and the power stations k 
were modeled as point masses ni A , m B , and mi, respectively. The power stations 
are connected with tether segments of lengths L^, as shown in Fig. 1. The tether 
segment Z* connects masses m*. and rrik+i. 


L« L, ... L v 


m A mi m 2 ... m k m k+ i m B 


Fig* 1. Structure of the momentum exchange tether system. 

All point masses and the masses of the tether segments arc assumed to be 
constant. Tether mass per unit length can vary along the tether. 

Positions of the tether elements with respect to a non-rotating geocentric 
reference frame OX Y Z are defined by the geocentric radius R as a function of the 
arclength .s measured along the unstretched tether from .4 to /?, and time t , 

R = R (s, t) y 

Positions of the end masses and embedded masses are 

R a = R(s > i,t) > Rb = R(sjj,<), Rfc = R 
The tension vector T of a perfectly flexible tether is tangent to the tether line, 

T = Tr, , " 7 = |R'I, ( 1 ) 

7 

where T is a unit vector tangent to the tether line, prime denotes differentiation 
with respect to the arclength s, and 7 is the local tether elongation. 

The tether tension T can be expressed as a function of the elongation 7, 
elongation rate 7, temperature 0, and other factors, 

T — T(s, t, 7, 7, < 9 ,...). ( 2 ) 
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k 


A 


Fig. 2. Positions of the tether system elements. 

Equations of motion of the tether system derived in [1] include ordinary and 
partial differential equations. The motion of the end masses and embedded masses 
is described by the ordinary differential equations 

m A = T 4 + m A g 4 4- F a 

m B R s = -T w 4- m B g B 4- F B (3) 

m k R-k = T fc+ - Tjt_ 4- m k g* 4- F* 

where dots denote differentiation with respect to time, g A , g B , and g fc are the 
gravity accelerations at points A, B, and k, respectively, while F 4 , F B , and F*. 
are non-gravitational forces acting on the end masses and embedded masses. The 
tether tension vectors are taken at the following points: at point A, T B at point 

B, T*_ at point k of segment I*_i, and T* + at point k of segment L k . 

The motion of the tether is described by the partial differential equation 

pR = T +pg + F (4) 

where dots denote differentiation with respect to time t , and primes denote differ- 
entiation with respect to the arclength s, p is the tether mass per unit length, and 
T is the tether tension. 

In the general case, the tether mass per unit length p can vary along the 
tether, but we will concentrate on a more practical case when p is constant along 
each tether segment connecting the end masses, 

P = Pk, k = l,...,K. (5) 
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2. EQUATIONS OF MOTION IN MINAKOV’S FORM 


Equations of motion in Minakov’s form, introduced into space tether dynamics 
in [2], are obtained as follows. First, we differentiate equation (4) with respect to 
the arclength s, 

pR' = T" +pg' + F'. 

Keeping in mind that R' = 7 x and T = Tx, we find that 


P (-yT + 27T + 7x) = T"x + 2 TV + Tt" + p g' + F'. 


( 6 ) 


By definition, the unit tangent vector x satisfies the condition 


The first differentiation with respect to the arclength and time yields two orthogo- 
nality conditions, 

(t,t) = 0, (t,x') = 0, (8) 

while the second differentiation yields two kinematic relations, 

(x, x) = -(x, x), (x, x") = -(x', x'). (9) 


Multiplying equation (6) by vector x (scalar product) and taking into account 
relations (7)-(9), we obtain a scalar equation of the second order with respect to 
tension, 

P7 - P7 (t, t) - T" - T(x',x') + p(g',x) + (F',x). (10) 

From this equation, we express the second derivative of tension as 


T " = +P1 ~ Pli*,*) -p(g',x) - (F',x). (11) 

Substituting expression (11) into equation (6), we arrive at the following equation 
of the second order with respect to the unit tangent vector X, 

p [ 7 T + 27 X + 7 t (t, x)] = T [x" + (x' , x 1 ) x] + 2 TV + pg' + F' - (pg' + F' , x) x. ( 12 ) 


This equation describes transverse motion of the tether and it has an obvious 
wave structure with a transverse wave velocity equal to 


Equation (10) describes longitudinal motion of the tether and it also has a 
wave structure. To reveal this structure, we consider a case of linear elasticity, 

T = T 0 + E( 7-7,), (14) 

where T 0 is a reference tension, E is the tether longitudinal stiffness, and 7, is 
he equilibrium elongation under the reference tension T 0 . According to (14) and 
keeping in mind that 7* may vary with time and temperature, we have 

- f „ 

7 “ E + 7 *' 

After substituting this relation into equation (10), we can clearly see that this 
equation has a wave structure with a longitudinal wave velocity equal to 


(15) 



The longitudinal wave velocity (15) is much higher than the transverse wave 
velocity (13), and therefore, the system of equations (10) and (12) is stifT in the sense 
that equation (12) describes high frequency longitudinal oscillations, while equation 
(10) describes transverse oscillations of much lower frequencies. The advantage of 
Minakov’s form of the equations lies in the separation of the high and low frequency 
wave structures in the tether dynamics equations. 

In the Newtonian field, 


the gravity gradient along the tether is calculated as 


g = - 




R 3 


r_Bi 

R 3 


+ ifHR>R') = ^[3e*(e*,T)- T ], 


(16) 


where fi B is the gravitational constant of the Earth, e K = R/R is a unit vector along 
the geocentric radius vector R, and R' = 7 T. The corresponding gravitational term 
in the equation of longitudinal motion (10) is reduced to 


(S\ r )= ~^f-l3(e R ,T) 2 -1], ( 17 ) 

while the respective gravitational term in the equation of transverse motion (12) 
takes the form 


S' ~ (s',t)t= ^ 3(e H ,T)[e R - (e*,T)T]. (18) 
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3. BOUNDARY CONDITIONS 


The equations of tether motion must be complemented with boundary condi- 
tions at the ends of each tether segment. To derive these boundary conditions, we 
express the accelerations of the end masses A and B and the embedded masses* m k 
from equation (3) 


n _ T* + F 4 

— _ b g4, 


Rb = 


m A 

— T g + F B 

m B 


+ gs, 


T) T *+ - T*_ -|- 
R ‘ = +g ‘> 


and from equation (4) at the ends of the tether segments, 


‘-Ft 1 ) 


+ g4, 


+ g B, 


+ g * 




+ gfc* 


This yields the following boundary conditions 


T f |a = Q,t, T'| b = Q b , 

T'| fe _ = Qt_, T'|*+ = Q k+> 


(19) 


where 


Q* = ~~ (T4 + F,*) - F| a , 

m A 

Qb = —(-T b +F b )-F\ b , 

TTlfi 

Q*- = ^ (T* + - T*_ + Ft) - F|*_, 

m t 

Q*+ = — (T* + - Tt_ + Ft) - F|*+. 

mt 


( 20 ) 


Note that the gravitational accelerations canceled out and are not included 
in the boundary conditions. 



Keeping in rnind that T' = T't + Tt\ according to (1), we multiply each of 
the relations {19} by the corresponding unit tangent vector t (scalar product) and 
arrive at the following conditions 

r'u = (Q A ,T < ), ru = (Q B , Ti ), 

r\ u . = (q».,t*-), tv - (Q 4 + ,T i+ ), (21) 

Relations (21) serve as boundary conditions for equation (10) describing longitudi- 
nal motion of the tether. 

Substituting relations (21) into (19), we find boundary conditions for equation 
(12) describing transverse motion of the tether, 


T '\a = ~ [Qa - (Qa,t a )t a ], 

j A 

T 'U = jT [Q« - 

T 'i*- = [Qfc- - (QA.,T fc _)r fc _], 


Ah+ = 


jrHQ*+ - (Qjb4,Tfc + )T fc+ ]. 


( 22 ) 


The boundary conditions are simplified when the non -gravitational forces F 
are neglected. The boundary conditions for the t ransverse motion of the tether are 
reduced to 

Aa = o, 

As = 0 , 

Ah- = ^ (T fe+ - **T fe _), (23) 

T 'l‘+ - P ir -*»-), 

Wfe -lk+ 

where 

= (Tfe+,Tjfc_), 


while the boundary conditions for the longitudinal motion of the tether take the 
form 


T\ a = — 1 

m A 

^ 




T'U = -^-T. 


m B 




m (*kTk+ - T fe _), T% + = 

mj. 


(24) 
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4. MOTION OF THE CENTER OF MASS 


The equations of tether motion derived in the previous sections describe evo- 
lution of the tether orientation and shape. They must be solved together with one 
of the ordinary differential equations (3) describing the orbital motion of one of the 
end masses or embedded masses, or with the equation of motion of the center of 
mass 

= (®o + 4>o), (25) 

where R c is the geocentric radius vector of the center of mass, M is the total mass 
of the tether system, G 0 is the sum of the gravitational forces 

rB 

G 0 = m A g„ + m B g B + ^ m/bgfc + / pg ds , (26) 

k 

and <|) 0 is the sum of the non-gravitational forces acting on the tether system, 

r B 

<I> 0 = F 4 F b + Y, F* + / F ds. (27) 

k 

As shown in the first part of this report [1], the sum of the gravitational forces 
can be represented as a series 


oo 

Go = Mg“ + (28) 

n=2 

where g° is the gravitational acceleration at the center of mass, g" are derivatives 
of the gravitational acceleration along the “mean” tether line R = R c + (.? - Sc ) Tl 
drawn through the center of mass, 


n= 1 

” c n! ds n 


and quantities I n represent the high order moments 

rB 

I n — fn A (sx — «c) n + (s B — Sc) n + ^ m* ( Sk — $c) n + / p(s — sc) n ds. 

t Ja 
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The direction of the “mean’ tether line T t is defined as follows. As shown in 
the first part of this report [1], the tether shape can be represented by the series 


R(.s,t) _ R c (t) + (s — «c)T|(f) + ^ q>t(0^n(s), (29) 

»= 2 

where q n (t) are generalized coordinates and U n (s) are the eigenforms of tether 
oscillations. Differentiating this equation with respect to the arclength, and keeping 
in mind that R/ = ^x, we derive 


T » =TT- J]q„(<)^n(*)- (30) 


The orthogonality conditions for the eigenforms can be expressed as 

J A PVlV' j ds = 0, t ^ j, (31) 

where P = P($) represents a normalized tension profile, defined as the solution to 
the following boundary problem [1], 


P ' = -/>(«- s c ), 

Pa = ~m A (s A - s c ), 

Pb = m B (s B - s c ), 

Pk+ = Pk- - m k (sk - s c ). 


(32) 


The first eigenform describes the rigid rotation £/, = s-s c , for which U' = 1. 
Multiplying equation (30) by PU[ = P , integrating over the entire tether length, 
and applying the orthogonality conditions (31), we find that 

B 

Pyrds. (33) 

i 

The integral of P , according to [1], is equal to the moment of inertia of the unde- 
formed tether system about the center of mass 



/ Pds = \\U l \\ = J c = 

J A 
f B 

m A (sa ~ s c ) 2 + m B (s B - Sc) 2 + ^ m* (s* — s c ) 2 + I p(s — s c ) 2 ds, 

k Ja 


(34) 
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and therefore, relation (33) can be rewritten as 


Ti 


Jc 



Pyrds. 


(35) 


Formula (35) defines the “mean” tether direction for any given tether shape under 
the assumptions of the current formulation. 

Knowing T, , we can calculate derivatives (26). As shown in [1], the Newtonian 
term of the gravitational field gives 


8c = 


PE 7," 

Rl +n 


[ e fl P n+l( x ) + e > ^n( a: )] 


(36) 


where c R = R C /|R C | is the unit vector of the geocentric direction to the center of 
mass, 7 ! = |Ti|, e x = Ti/h, P n {x) are Legendre polynomials, and x = ~(e fi ,ei). 

To close the system of equations, we need to find the geocentric positions of 
all elements of the tether system. For a given tether shape, we know the relative 
positions r = R — R* of all elements, and from the definition of the center of mass, 
we can find its relative position 

r c — + 'y + I prdsj. (37) 


Using this information, we can calculate the geocentric position of the end 
point R a — R c — t c and the geocentric positions of all other points of the tether 
system R = R^ + r. 


5. QUASI-STATIC TENSION 


The longitudinal oscillations of the tether system are directly affected by the 
internal friction in the tether, while the transverse oscillations are are affected 
indirectly, through a weak coupling with the longitudinal oscillations. 

As discussed in [2], one way to model internal friction in the tether is to add 
a term px 7, where \ is an effective damping coefficient, to the left side of the equa- 
tion of longitudinal motion (10). It is shown in [2] that the energy dissipation in 
the tether can rapidly damp out high frequency longitudinal oscillations, and after 
a short period of t im e the longitudinal motion becomes quasi-static. Tension vari- 
ations in the quasi-static motion are induced by relatively slow transverse motions 
and can be described by equation (11) without the term p 7 , 

r - T(V, t') - - p (g'> t) - (F». ( 38 ) 
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At any given moment f, this equation can be treated as an ordinary differential 
equation with respect to the tether tension T. Along with the boundary conditions 
i 21 ?, it defines the tether tension profile as a function of the dynamic state of the 

tether, including the tether shape, elongation and the angular rate distribution 
along the tether. 

Similar to (24), the boundary conditions ( 21 ) can be presented as 


T'U = ^T i + / A , 

m A 

T'\ B = -£*- Ta +f B , 

Til 

T> ^+ = ( Tk + - + fk+, 

Tflfc 


(39) 


where 


ffl A 

* = r.,T»)-(r | „T,). 

/*-= ~ (F»,t,-)-(FU-,t*_), 

At = (Fi>1»+) -(F|t+,T, + ). 


As we see in relation ( 16 ) , the gravity gradient g' calculated along the tether 
is proportional to the elongation 7 , g' - 7 g'. This is true for any gravitational 
field, not only Newtonian. Also, it has been shown in [ 2 ] that the aerodynamic 
and Ampere forces, as well as solar radiation pressure acting on a tether element 
are proportional to the elongation, F = 7 F 0 . Formally, the gradient of these forces 
F' will include quadratic terms in 7 because of the variation of the environmental 
parameters along the tether, however, these terms can be linearized with respect 
to small variations of elongation. 

Under these assumptions, if we use the linear elasticity relation (14) between 
the elongation and tension, then we will have 


7 = 7 , + (T - T 0 )/E, (40) 

and equation (38) will be linear with respect to T . The boundary conditions ( 39 ) 
will also be linear with respect to T. 

The general solution of the linear boundary problem (38)— f 39) at any given 
moment t can be represented as a sum of two solutions, 

^ (-SO — ^o( s jt) + CiTi(s,t). (41) 
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Solution T 0 is obtained as follows. We set T 0A to zero and determine T' from 
the first equation of the system (39), 

T 0 A = 0, Tq A = f A . (42) 

1 hen, we integrate equation (38) with the initial conditions (42) over the first tether 
segment, and find and T' k _ at the end of the segment. Now, we nse the third 

equation of the system (39) to find the tension at the beginning of the next tether 
segment 

Ti+ = 7T k [ r ‘- + jt (r *- “*-)]• («) 

We then substitute quantity (43) into the last equation of the system (39) to find 
the derivative at the beginning of the next segment. 

This process is repeated until we arrive at the end B and determine the values 
of T cb and T‘ IB . Generally, these values will not satisfy the second equation of the 
system (39). We need solution T[ to match the boundary condition at the end B. 

Solution 7\ is calculated as follows. We introduce a new system of equations 

(38’ j (39') which retains only terms linear in T from the original equations (38)- 

(39). We set T lA to some reference tension , 

T,a=T„ ( 44 ) 

find T; a from the first equation of the homogeneous system (39'), and integrate the 
homogeneous equation (38') over the first tether segment. We then use the third 
and the last equation of system (39') to find 2*+ and T‘ k+ at the beginning of the 
next segment, and continue integration until we reach the end B. 

Now, we can substitute the general form of solution (41) into the second 
equation ol (39) and find the coefficient Cj from the following equation 

+ m. *B Tr3B + Cl ( T ' b + ^ TiB ) = /b ’ ( 4S ) 

One of the major benefits of this approach is that we reduced the initially stiff 
system of equations (10), (12) to a non-stiff system (12), (38) by abstracting from 
the fast transient processes in longitudinal oscillations. This does not mean that 
we neglected tether elongation, we only filtered out its high frequency components. 


6. SIMULATION APPROACH 


For a numerical solution, equations (12), (22), (38), (39) can be discretized 
with respect to the arclength s. This will reduce the partial differential equation 
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(12) to a set of ordinary differential equations at the discretization nodes. These 
equations are complemented with the ordinary differential equation (25) describing 
the motion of the center of mass. 

The number of discretization nodes is determined by the required accuracy. 
To get an idea of a typical number of nodes, let us consider a very simple model of 
a siring clamped at the ends A and H. Its transverse oscillations are decribed by 
the string equation 

pil — Tu'\ 

where p is the string mass per unit length, T is the tension, assumed constant 
along the string, u = u(s,f) is the transverse displacement, dots denote differenti- 
ation with respect to time i, and primes denote differentiation with respect to the 
ardength s. 

The boundary conditions, are, obviously, 

Ua — Ub = 0 , 


A discretized string equation looks as follows 

2 u k -\- 1 ~ + Ufa— i 


u = v: 


as* 


where v t - y/Tjp is the transverse wave velocity, ti* is the transverse displacement 
at the node fc, as is the distance between the neighboring nodes, 


AS = 


N * 


N is the number of nodes, and L is the length of the string. 

Partial solutions representing the natural modes of oscillations are found in 
the form 

u k - sin(J?„f) sin(Tm^), n = 1,2, . . . 

where O n are the eigenfrcquencies, and — k /N, k = 0, 1, , N . 

Substituting this form of solution into the discretized string equation, we find 
that the eigenfrcquencies of the discretized problem are equal to 

n ojir Vt • f‘ Kn \ v * (-, 

L \2N J L y 24 N 2 J 

The exact eigenefrcquency obtained from the original partial differential equa- 
tion is equal to 

rt* Vt 

il n — irn — * 

L 
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The relative difference is 


^n. ~ ^ 1 /irn\ z 

fl* n ~ 24 \JV"/ * 

Because of the difference in eigene frequencies, the phase difference between 
the simulated and exact solutions will accumulate with time as ay? = (Q n - £2' n )t. 
If we need to simulate the excitation of the n-th mode with a relative accuracy e T1 

over a period of time af, then the difference between the eigenefrequencies must be 
limited as 

IA. - All < % 

At 

To satisfy this relation, the number of nodes must be at least 


N > 

nj 



For a typical momentum exchange system with L = 90 km and v t = 1 krn/s, 
it takes at least 125 nodes to maintain an accuracy of 1% for the first mode of 
transverse oscillations over one 3-hour orbit. 

If an explicit integration scheme is used, then the time step will be limited by 
the condition of stability of the numerical solution. In general, this restriction has 
a form of 


M < kt min ( — — 1, 
* 1 ®ti J ' 


where k t is a numeric coefficient depending on the implementation, as; are the 
distances between the nodes, and vn are the maximum transverse wave velocities 
(13) between the corresponding nodes. 

By comparison, in a lumped mass model, the time step is limited by 


A f> ^ 1?^ 



where k r is a numeric coefficient depending on the implementation, and v ei are the 
maximum longitudinal wave velocities (13) between the corresponding nodes. 

The longitudinal wave velocity is much higher than the transverse wave ve- 
locity, and therefore the time step will be much smaller. This is why simulations 
based on the Minakov’s form of the equations of motion are much faster than those 
based on lumped mass models. 

Two numerical models were implemented, one using Minakov’s form of the 
equations, and the other using a lumped mass model. With a sufficient number of 
nodes, as explained above, the results of the simulations were consistent between 
the models within the required accuracy. As expected, the Minakov’s type model 
was by an order of magnitude faster than the lumped mass model. 
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7. RESONANT EXCITATION 


It has been noted in simulations that the transverse oscillations of a typical 
momentum exchange system go through cycles of a long-term amplitude modula- 
tion, when the amplitude first builds up during a number of consecutive perigee 
passages, and then decreases again. 

Fig. 3 shows a history of the transverse amplitude variations over the course 
of twelve orbits. The amplitude profile is shown in light gray, while the perigee 
passages are marked with black clusters underneath. 



Fig. 3. Resonant excitation of the transverse oscillations. 

We see that even without the initial excitation the amplitude of transverse 
oscillations steadily grows during five perigee passages until it reaches a relatively 
high level of over 120 m T stays at this level for another perigee passage* and then 
steadily decreases during the next five perigee passages. This cycle repeats with 
minor variations. 

An in-depth analysis of this phenomenon based on the modal formulation de- 
veloped in the first part of this report reveals that the mechanism of this excitation 
is resonant in nature and is caused by a sharp tuning between the natural frequency 
of the main mode of transverse oscillations U 2 and the frequency of the gravity gra- 
dient variation, which cycles twice per every rotation of the tether system about 
its center of mass. 

This timing is persistent in the sense that the parameters of the tether system 
must be changed radically to get away from the resonance. In particular* changing 
the spin rate has very little effect because the transverse frequencies are proportional 
to the spin rate, as shown in [1]. 
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The modal formulation of Part f accurately describes the motion of the tether 
system in general and the resonant excitation of transverse oscillations in particular 
when the term Q* in the main modal equation (37) derived in Part I is used in its 
general form (38). The simplifications of the term Q n demonstrated for a special 
(.ase o an ideally tapered tether should remain restricted to this special case and 

should not be applied in case when the MXER tether is not carrying the full load 
lor which it was designed* 

tor practical purposes, the modal solution requires very few modes and is 
computed much faster than the numeric solutions based on the Minakov’s formu- 
latiGii, not to mention lumped mass models. 


8* SENSITIVITY TO NON-GRAVITATIONAL PERTURBATIONS 
8.1* Method of Computation. 

To evaluate the sensitivity of the motion of the momentum exchange system to non- 

gravitational perturbations, we use the modal decomposition method described in 
Part I of this study [1]. 

As shown in [1], the effect of non-gravitational forces on the excitation of the 
eigen form V n is determined by a generalized force 

cp n = F A U nA + F B U nB + U nk + f FU n ds. (46) 

k Ja 

It is computationally expensive to calculate the integral of F U n over the entire 
tether length on each step. To make this calculation more efficient, we will use the 
following quadratic approximation on each tether segment s k < s < s k+1 

= (2f -e)F + (1 - 4£ 2 )F(s fe+x/2 ,f) + (f + 2f a )F (s k+u t), (47) 

where £ (s f 2 )/(■**-(- 1 — an d Sfe-1-1/2 = (•*& T a *+i)/2 is the middle of 

the segment, so that -1/2 < £ < 1/2, 

Using approximation (47), we have 


f FUnds = y 

JA t 


yZ Qnk, 


where 

rk+l 

Qnfc=/ F U n ds = 

J k 

i 2 *nk - 4iiJF(sfc,f) + {I° k - 4Il k )F(s h+1/2> t) + (/^ + 2/^ A .)F(s fc+ i,t), 


53 



and 


T° 

*nk 


rk-n 

Jk 


U n ds, 


rM-i 

= L 


f- J ' n f 


1 2 


nk 


fk+l 

Jk 


U n ( 2 ds. 


SSSAZZZZZ' - ~ d — - =s 

With tjpical parameters, the accuracy of this presentation of the generalized 

forces meets ° r as «**«-«*. ** *££££ 
dynamics of the momentum exchange system. 

8.2. Aerodynamic Forces. 

matedT" b 121 ’ tl,e “ rodynMr,ic for " “‘“g on » ‘other segment can be appro*!- 


P = -PtldtJ 


|T xv.|( (l + |) v a - 1 £Vt j +*/(!- e)Jv w j(v a - v r ) , 


(49) 


where d t is the tether diameter, p a is the air density, v c is the tether velocity relative 

a °, hc "5 Vr " ( v - i y) T 18 the component of the velocity v ft along the tether line 
and e and i/ are small parameters, 0 < e, u < 0 1 

The aerodynamic forces acting on the end masses and the embedded masses 
Ere calculated as for regular satellites [4]* 

, |’ or ^ lloytether [5] consisting of separate strands kept apart from each other 
formula (49) must be applied to all individual strands. One must also take into 
account the overshadowing of the strands. In general, a Hoytethcr will have a larger 
total exposed area of the tether and a higher air drag. 

With typical system parameters and Hoytethers having between 16 and 24 
strands, it has been determined in numerical simulations that a 5% variation in the 
aerodynamics forces due to air density variations or uncertainties of the aerody- 
namic parameters of the tether will cause a 2 m shift of the tether tip A after one 
orbit . 


8.3. Ampere Forces. 

Tlie Ampere force acting on a tether segment can be represented 

F = Ijt x B 


as 


(50) 


vector 1 ' 5 th<! elCC ‘ riC " lrient “ ‘ hC tether ’ "" d B “ the ^““gnetic induction 
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Normally, Ampere forces should not be applied during a preparation for a 
rendezvous. However, to evaluate the sensitivity to the variations of the Ampere 
forces, we will consider a purely hypothetical situation when there is an electrical 
“leak” in the system, which results in a very small alternate electric current of a 
1 mA amplitude flowing in the conductive segments of the tether in the direction 
of the EMF and modulated proportionally to the EMF. It was also assumed for 
simplicity that the “leak” current shuts off at altitudes above 2000 km for the lack 
of electrons. 

It has been observed in simulations that this electric “leak” during the perigee 
passage causes a 3 m shift of the tether tip A after one orbit. 


8.4. Solar Radiation Pressure. 


The force of the solar radiation pressure acting on a round tether can be approxi- 
mated as [2] 


F = -p t d t 7 |t x e,| 


V, 

4 . , 1 

i( i+ j) 

e, - - *(e,,T)T 


(51) 


where e, is a unit vector of the direction to the Sun, p 3 is the solar radiation 
pressure, and x is the reflection factor. 

The solar radiation pressure forces acting on the end masses and the embedded 
masses are calculated as for regular satellites [4] . 

For a Hoytether [5] consisting of separate strands kept apart from each other, 
formula (51) must be applied to all individual strands, taking into account their 
overshadowing. In general, a Hoytether will have a larger total exposed area of the 
tether. 

With typical system parameters and Hoy tethers having between 16 and 24 
strands, it has been determined in numerical simulations that a 5% variation in the 
forces of solar radiation pressure due to uncertainties of the reflective and geometric 
parameters of the tether will cause a 1 to 4 m shift of the tether tip A after one 
orbit. 


8.5. Thermal Expansion. 

Thermal expansion changes the length of the tether segments and thus changes 
the moment of inertia of the momentum exchange system. The variations of the 
moment of inertia result in variations of the tether spin rate. 

The moment of inertia is proportional to the second power of the average 
elongation 7, , and in a simple case of a free rotation, when the angular momentum 
is conserved, the angular velocity Q satisfies the relation 

7?^ = 7o^°* (52) 
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When the average elongation 7, varies with the temperature 0 as 


7i=7o + ot t a 0 , (53 ) 

where a t is the thermal expansion coefficient, then the angular rate will vary as 


7o 


7o 


If the average deviation of the actual temperature from the modeled temper- 
ature is a<9 ot , during a time af, then the angle of the tether in-plane orientation 
will deviate by ai? = and the transverse displacement of the tether end ,4 

will be estimated as 

„ r 2a ( a(9 at , 

Al) L AC fti — - do At L ac , 

where L AC is the distance from the end A to the center of mass C. 

According to [3], Zylon, a candidate tether material, has a negative thermal 
expansion coefficient of a, = -6x 10“ 6 1/K. With a typical spin period of 6.3 min 
and a distance from the end to the mass center of 75 km, we find that after one 
orbital period of M =: 3 hours, the tether end will drift away from its anticipated 
position by about 163 m per each Kelvin of the average temperature deviation. 
This sensitivity is attributed to the high spin rate Q 0 . 

8 .6. Creep. 

The average elongation 71 will gradually increase with time due to creep, 

7i=7o + « e ^, (54) 

where a c is the creep coefficient. In the approximation (52), the average spin rate 
will be dropping as 

/o 2^7, 2 a c At „ 

« - — — o Q = — — n 0 . 

7o 7o 

and the deviation of the angle of the tether in-plane orientation will be growing as 

ad fe - — — 

7o 

The transverse displacement of the tether end A will build up with time as 

< 1 ^ L ac ^ L A c } 

7o 

where l A c is the distance from the end A to the center of mass C. 
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As reported in [3], the creep coefficient a c of Zylon-HM under high stress 
is on the order of 10^ fi 1/hour, For a 100- km tether, this will result in tether 
elongation at a rate of 2*4 m/day* With a typical spin period of 6.3 min and a 
distance from the end to the mass center of 75 km, we find that after one orbital 
period of aI = 3 hours, the tether end will drift away from its anticipated position 
by about 41m. For Zylon-AS, the creep coefficient and the tether end displacement 
are 3 times higher* As with thermal expansion, this sensitivity is attributed to the 
high spin rate 

8,7* Mass Loss. 

There are several processes that will contribute to the variation of the mass of the 
momentum exchange system: 

(1) outgassing 

(2) sublimation 

(3) micrometeorite damage 

(4) hollow cathode emission 

(5) molecular deposition (mass addition) 

While there is no reliable data on most of them, it has been estimated that a typical 
momentum exchange system may lose between 5 and 150 grams of its mass per day. 
As a result, of this process, the center of mass of the tether system will be slowly 
shifting along the tether. 

Let us consider the motion of a system with a straight, quasi-rigid tether, 
whose elements m* are positioned at points along the line 


Ri = Re? + (■Si — (55) 

where x is a unit vector, and subscript C refers to the center of mass. By the 
definition of the center of mass, 

5^ mi ( si - s c ) = 0. (56) 

i 

According to (56), when the masses of the elements change, the center of mass 
slides along the tether at a rate of 

*c = m (57) 

i 

where M is the total mass of the tether system, 

M = ^ rrti. 

I 
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Differentiating (55) with respect to time, we obtain 


R, — Re — &c t — 2 Sc t •+• (it — ic) t. (58) 

Multiplying (58) by m i and performing summation, we find that 

^miRi = M(R C - ict — 2i c t), (59) 

ft 

where the terms with T canceled out because of the definition of the center of mass 
(56). 

For each element of the system, 

m»Ri = T, + F, + W it (60) 

where T,, F,-, and W,- are the sums of all internal, external, and reactive forces 
acting on this element. Substituting equations (60) for the individual elements into 
(59), we arrive at the following equation of motion of the center of mass 

R-c = ic T + 2 i c T + — ^^(F< + Wj). (61) 

ft 

All internal forces T* cancel out after summation, as dictated by the Newton^s third 
law. The terms with the derivatives of s c describe the perturbation of the orbital 
motion because of the relocation of the center of mass* 

Multiplying (58) by T (vector product) and performing summation, we find 

that 

Jc ft = 5>i -«c)tx(F, + W ( ), (62) 

ft 

where J c is the moment of inertia about the center of mass 

J c = - sc) 2 , 

t 

and O is the angular velocity of the tether rotation 

Q = T X T. 

It is interesting to note that even though the moment of inertia J c is changing 
because of the mass loss, it does not directly result in the variation of the angular 
velocity, unlike the case of creep and thermal expansion. 

While the hollow cathodes are expected to create a small torque of reactive 
forces Wj, it is not obvious that the other processes involved in the mass exchange 
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with the environment will be anisotropic enough to create any noticeable torque of 
the reactive forces. 

Assuming that the reactive torques are negligible, the gravitational field is 
Newtonian, and the non-gravitational forces are small, the in-plane perturbation of 
the orbit of the center of mass can be described, according to (61), as 


Sx — 2u) Sy — tb Sy — 


Sy -I- 2u ; Sx + cj Sx — 



2 /*g \ 

R%) 

R%) 


Sx = s c cos 0 — 2 s c f? sin t?, 
Sy = $c sin i? 4* 2sc cos i?, 


(63) 


where Sx and Sy are the vertical and horizontal displacements of the, center of mass, 
u) is the orbital angular rate, f? is the spin rate, and t? is the angle between the 
tether and the local vertical. For a fast spinning tether, J? u;, we have i? % f2t, 
and the order of magnitude of the perturbation is estimated as 


Sx 2 + Sy 2 



These displacements are expected to be very small (less than a millimeter) because 
of the fast spin and slow shift of the center of mass. 


9. ESTIMATION AND CONTROL REQUIREMENTS 


As shown in the previous section, the position of the tether tip in the momen- 
tum exchange system is very sensitive to variations of the system and environmental 
parameters, when considered from the standpoint of the rendezvous precision re- 
quirements. 

The variations producing the most impact on the tether tip positioning are 
related to the lengths of the tether segments, as we have seen in the case of thermal 
expansion and creep. 

If the distances between the neighboring tethered modules could be measured 
with an accuracy of a few millimeters, then we could use this information to sig- 
nificantly improve the estimate of the current state of the momentum exchange 
system. 

Every dynamic process in the momentum exchange system, one way or an- 
other, leaves its signature in the longitudinal motion of the tether. Using a detailed 
dynamic model, we may be able translate a continuous stream of precise distance 
measurements, combined with other data, into a reasonably accurate estimation of 
the tether system state. 


Furthermore, using the information about the system state, we could com- 
pensate for some unexpected (or unmodeled) deviations from the projected path. 
In particular, deviations of the rotation about the center of mass could be corrected 
by slightly varying the tether length(s) or changing the mass distribution of the end 
bodies or power stations. 

Given extremely short rendezvous windows with very fast relative motions, 
and various inherent uncertainties of the environment, it is very clear that elaborate 
estimation and control algorithms must be developed to make successful rendezvous 
with the momentum exchange system possible. 


10. CONCLUSIONS 


The task of prediction of the motion of a momentum exchange system with 
a 9-digit accuracy is extremely challenging, and it is relentlessly testing our ability 
to gain fundamental insights into the nature of the tether dynamics in this system. 

It has been shown that the modal decomposition approach developed in the 
first part of this study is a very powerful and precise tool for the simulation of the 
dynamics of momentum exchange tethers. 

However, to get a prediction with the required accuracy using this tool, one 
must have precise inputs, including the initial state, system parameters, and the 
environmental models. In simulations, the tether tip positioning has been observed 
to be quite sensitive to even small variations of the parameters involved in the 
calculations. 

It is therefore imperative that precise estimation and control algorithms be 
developed to compensate for the inevitable uncertainties and support high precision 
rendezvous. 
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